── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks plotly::filter(), stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)library(tidyr)library(zoo)
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Importing Data
## data extracted from New York Times state-level data from NYT Github repository# https://github.com/nytimes/covid-19-data## state-level population information from us_census_data available on GitHub repository:# https://github.com/COVID19Tracking/associated-data/tree/master/us_census_data### FINISH THE CODE HERE #### load COVID state-level data from NYTif (!file.exists("us_census_2018_population_estimates_states.csv"))download.file(url ="https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv",destfile ="us_census_2018_population_estimates_states.csv",method ="libcurl",timeout =60 )state_pops <- data.table::fread("us_census_2018_population_estimates_states.csv")
if (!file.exists("us_states.csv"))download.file(url ="https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv",destfile ="us_states.csv",method ="libcurl",timeout =60 )cv_states <- data.table::fread("us_states.csv")
We need to change the variable type of states and abb. We also need to make the date an actual date we can pull from rather than just an IDate.
# format the datecv_states$date <-as.Date(cv_states$date, format="%Y-%m-%d")# format the state and state abbreviation (abb) variablesstate_list <-unique(cv_states$state)cv_states$state <-factor(cv_states$state, levels = state_list)abb_list <-unique(cv_states$abb)cv_states$abb <-factor(cv_states$abb, levels = abb_list)### FINISH THE CODE HERE # order the data first by state, second by datecv_states = cv_states[order(cv_states$state, cv_states$date),]# Confirm the variables are now correctly formattedstr(cv_states)
# Inspect the range values for each variable. What is the date range? The range of cases and deaths?head(cv_states)
state date fips cases deaths geo_id population pop_density abb
1: Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
2: Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
3: Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
4: Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
5: Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
6: Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
summary(cv_states)
state date fips cases
Washington : 1158 Min. :2020-01-21 Min. : 1.00 Min. : 1
Illinois : 1155 1st Qu.:2020-12-06 1st Qu.:16.00 1st Qu.: 112125
California : 1154 Median :2021-09-11 Median :29.00 Median : 418120
Arizona : 1153 Mean :2021-09-10 Mean :29.78 Mean : 947941
Massachusetts: 1147 3rd Qu.:2022-06-17 3rd Qu.:44.00 3rd Qu.: 1134318
Wisconsin : 1143 Max. :2023-03-23 Max. :72.00 Max. :12169158
(Other) :51184
deaths geo_id population pop_density
Min. : 0 Min. : 1.00 Min. : 577737 Min. : 1.292
1st Qu.: 1598 1st Qu.:16.00 1st Qu.: 1805832 1st Qu.: 43.659
Median : 5901 Median :29.00 Median : 4468402 Median : 107.860
Mean : 12553 Mean :29.78 Mean : 6397965 Mean : 423.031
3rd Qu.: 15952 3rd Qu.:44.00 3rd Qu.: 7535591 3rd Qu.: 229.511
Max. :104277 Max. :72.00 Max. :39557045 Max. :11490.120
NA's :1106
abb
WA : 1158
IL : 1155
CA : 1154
AZ : 1153
MA : 1147
WI : 1143
(Other):51184
min(cv_states$date)
[1] "2020-01-21"
max(cv_states$date)
[1] "2023-03-23"
The date range is from Jan 21, 2020 to March 23, 2023. Cases range from 1 to 12169158 and deaths range from 0 - 104277.
Cases, Deaths, and Outliers
for (i in1:length(state_list)) { cv_subset =subset(cv_states, state == state_list[i]) cv_subset = cv_subset[order(cv_subset$date),]# add starting level for new cases and deaths cv_subset$new_cases = cv_subset$cases[1] cv_subset$new_deaths = cv_subset$deaths[1]### FINISH THE CODE HEREfor (j in2:nrow(cv_subset)) { cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j -1] cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j -1] }# include in main dataset cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths}
# Focus on recent datescv_states <- cv_states %>% dplyr::filter(date >="2021-06-01")
### FINISH THE CODE HERE# Inspect outliers in new_cases using plotlyp1<-ggplot(cv_states, aes(x = date, y = new_cases, color = state)) +geom_point(size = .5, alpha =0.5)ggplotly(p1)
#p1<-NULL # to clear from workspace
California, Texas and Florida have high outliers. Florida, Colorado, Pennsylvania, Tennessee, Kentucky, Texas, Washington, Nebraska, and other states have outliers because they have data points for new cases that are below zero which is not possible.
p2<-ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +geom_point(size = .5, alpha =0.5)ggplotly(p2)
#p2<-NULL # to clear from workspace
# set negative new case or death counts to 0cv_states$new_cases[cv_states$new_cases<0] =0cv_states$new_deaths[cv_states$new_deaths<0] =0
# Recalculate `cases` and `deaths` as cumulative sum of updated `new_cases` and `new_deaths`for (i in1:length(state_list)) { cv_subset =subset(cv_states, state == state_list[i])# add starting level for new cases and deaths cv_subset$cases = cv_subset$cases[1] cv_subset$deaths = cv_subset$deaths[1]### FINISH CODE HEREfor (j in2:nrow(cv_subset)) { cv_subset$cases[j] = cv_subset$new_cases[j] + cv_subset$cases[j -1] cv_subset$deaths[j] = cv_subset$new_deaths[j] + cv_subset$cases[j -1] }# include in main dataset cv_states$cases[cv_states$state==state_list[i]] = cv_subset$cases cv_states$deaths[cv_states$state==state_list[i]] = cv_subset$deaths}# Smooth new countscv_states$new_cases = zoo::rollmean(cv_states$new_cases, k=7, fill=NA, align='right') %>%round(digits =0)cv_states$new_deaths = zoo::rollmean(cv_states$new_deaths, k=7, fill=NA, align='right') %>%round(digits =0)
# Inspect data again interactivelyp2<-ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +geom_line() +geom_point(size = .5, alpha =0.5)ggplotly(p2)
#p2=NULLtable(is.na(cv_states))
FALSE TRUE
377419 673
Step 5
### FINISH CODE HERE# add population normalized (by 100,000) counts for each variablecv_states$per100k =as.numeric(format(round(cv_states$cases/(cv_states$population/100000),1),nsmall=1))cv_states$newper100k =as.numeric(format(round(cv_states$new_cases/(cv_states$population/100000),1),nsmall=1))
# add a naive_CFR variable = deaths / casescv_states = cv_states %>%mutate(naive_CFR =round((deaths*100/cases),2))# create a `cv_states_today` variablecv_states_today =subset(cv_states, date==max(cv_states$date))
Scatterplots
# pop_density vs. casescv_states_today %>%plot_ly(x =~pop_density, y =~cases, type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
Warning: Ignoring 1 observations
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Pop Density vs Cases after Filtering
# filter out "District of Columbia"cv_states_today_filter <- cv_states_today %>%filter(state!="District of Columbia")# pop_density vs. cases after filteringcv_states_today_filter %>%plot_ly(x =~pop_density, y =~cases, type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
Warning: Ignoring 1 observations
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Pop Density vs Deaths per 100k
# pop_density vs. deathsper100kcv_states_today_filter %>%plot_ly(x =~pop_density, y =~deathsper100k,type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
Warning: Ignoring 1 observations
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Adding Hoverinfo
# Adding hoverinfocv_states_today_filter %>%plot_ly(x =~pop_density, y =~deathsper100k,type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5),hoverinfo ='text',text =~paste( paste(state, ":", sep=""), paste(" Cases per 100k: ", per100k, sep="") , paste(" Deaths per 100k: ", deathsper100k, sep=""), sep ="<br>")) %>%layout(title ="Population-normalized COVID-19 deaths (per 100k) vs. population density for US states",yaxis =list(title ="Deaths per 100k"), xaxis =list(title ="Population Density"),hovermode ="compare")
Warning: Ignoring 1 observations
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Step 7: Explore trend interactively
### FINISH CODE HEREp <-ggplot(cv_states_today, aes(x=pop_density, y=deathsper100k, size=population)) +geom_point() +geom_smooth() +labs(title ="Scatterplot of pop_density vs. newdeathsper100k", x ="Population Density",y ="New Deaths per 100k")ggplotly(p)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: size
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
I don’t think population density is a correlate of new deaths per 100k.
Step 8: Multiple Line Charts
### FINISH CODE HERE# Line chart for naive_CFR for all states over time using `plot_ly()`plot_ly(cv_states, x =~date, y =~naive_CFR, color =~state, type ='scatter', mode ='lines')
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
States that had an increase in September decreased over time and had other spikes in December.
### FINISH CODE HERE# Line chart for Florida showing new_cases and new_deaths togethercv_states %>%filter(state=="Florida") %>%plot_ly(x =~date, y =~new_cases, type ="scatter", mode ="lines") %>%add_trace(x =~date, y =~new_deaths, type ="scatter", mode ="lines")
There is around a month delay between peaks of cases and peaks of deaths.
The states that now stand out are Rhode Island, New York, New Jersey, and Alaska.
# Create a second heatmap after filtering to only include dates every other weekfilter_dates <-seq(as.Date("2021-06-15"), as.Date("2021-11-01"), by=14)
# Create a heatmap using plot_ly()plot_ly(x=colnames(cv_states_mat2), y=rownames(cv_states_mat2),z=~cv_states_mat2,type="heatmap",showscale=T)
Step 10: Map
### For specified datepick.date ="2021-10-15"# Extract the data for each state by its abbreviationcv_per100 <- cv_states %>%filter(date==pick.date) %>%select(state, abb, newper100k, cases, deaths) # select datacv_per100$state_name <- cv_per100$statecv_per100$state <- cv_per100$abbcv_per100$abb <-NULL# Create hover textcv_per100$hover <-with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))# Set up mapping detailsset_map_details <-list(scope ='usa',projection =list(type ='albers usa'),showlakes =TRUE,lakecolor =toRGB('white'))# Make sure both maps are on the same color scaleshadeLimit <-125# Create the mapfig <-plot_geo(cv_per100, locationmode ='USA-states') %>%add_trace(z =~newper100k, text =~hover, locations =~state,color =~newper100k, colors ='Purples' )fig <- fig %>%colorbar(title =paste0("Cases per 100k: ", pick.date), limits =c(0,shadeLimit))fig <- fig %>%layout(title =paste('Cases per 100k by State as of ', pick.date, '<br>(Hover for value)'),geo = set_map_details )fig_pick.date <- fig################ Map for today's date# Extract the data for each state by its abbreviationcv_per100 <- cv_states_today %>%select(state, abb, newper100k, cases, deaths) # select datacv_per100$state_name <- cv_per100$statecv_per100$state <- cv_per100$abbcv_per100$abb <-NULL# Create hover textcv_per100$hover <-with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))# Set up mapping detailsset_map_details <-list(scope ='usa',projection =list(type ='albers usa'),showlakes =TRUE,lakecolor =toRGB('white'))# Create the mapfig <-plot_geo(cv_per100, locationmode ='USA-states') %>%add_trace(z =~newper100k, text =~hover, locations =~state,color =~newper100k, colors ='Purples' )fig <- fig %>%colorbar(title =paste0("Cases per 100k: ", Sys.Date()), limits =c(0,shadeLimit))fig <- fig %>%layout(title =paste('Cases per 100k by State as of', Sys.Date(), '<br>(Hover for value)'),geo = set_map_details )fig_Today <- fig### Plot together subplot(fig_pick.date, fig_Today, nrows =2, margin = .05)
In 2021 there were a lot more states with many cases per 100k that were above 20 and even some that were close to 100. In 2023, most states have less than 10 cases of covid per 100k.